#------------------------------------------------------------------------------
# Import libraries
#------------------------------------------------------------------------------

rm(list = ls())
cat("\014")
library(tidyverse)
library(PanelMatch)
library(data.table)
library("zoo")
library("tictoc")
library('fst')
library("fixest")
library("PanelMatch") 
library("patchwork")
library("rio")
library("magrittr") 
library("janitor")
library('did')
library("panelView")
library("vtable")
library("RPushbullet")

#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
# Some comments
# This codes takes too much time so it should be run in the cluster. We bring 
# to you the final object after estimation to get figures very quickly. If you 
# would like to follow all the process you should run 
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
# Define Paths
#------------------------------------------------------------------------------

# R studio
setwd( dirname(rstudioapi::getActiveDocumentContext()$path) )
# R default : unccoment if you use default R
# setwd(getSrcDirectory(function(){})[1])
#------------------------------------------------------------------------------





#------------------------------------------------------------------------------
# Import data
#------------------------------------------------------------------------------

load( "panelmatch_ps_vcf_forest_index.RData" )
load( "fig7_stateMining_exact_forest_index_obj.RData" )
load( "fig5_stateMining_exact_forest_index_obj.RData" )

#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

vcf <- fread("vcf_data_complete.csv", sep = ",")
setnames( vcf, "d", "D")
vcf_data <- copy( vcf )

#------------------------------------------------------------------------------


#------------------------------------------------------------------------------
# Load data
#------------------------------------------------------------------------------

vcf = vcf_data[year >= 1995]
vcf[, t := year - 1995]
ex_ante_med = quantile(vcf$cover_1990, 0.5)
above_med = vcf[cover_1990 > ex_ante_med]
above_med[ , never_treated := max(D) == 0, cellid ]

controls_pre1 = above_med[never_treated == 1 & year < first_pesa_exposure, 
                          mean(forest_index)]
treat_pre1    = above_med[never_treated == 0 & year < first_pesa_exposure, 
                          mean(forest_index)]

#------------------------------------------------------------------------------

m0_ps <- feols( forest_index ~ 1| cellid + year , data = above_med )
m1_cbps <- feols( forest_index ~ 1| cellid + year , data = above_med )


#------------------------------------------------------------------------------
# Generate Tables
#------------------------------------------------------------------------------

# Change name of coefficients
treatmap =c(
  # GFC
  "def_ha"       = "Annual Deforestation in Hectares",
  "village"      = "Village",
  "village[t]"   = "Village + Village TT",
  # VCF
  "forest_index" = "Forest cover index",
  "green_index"  = "Non-Forest Vegetation",
  "built_index"  = "Bare Ground Indices",
  "cellid"       = "Pixel",
  "cellid[t]"    = "pixel + pixel TT",
  # both
  "yr"           = "Year",
  "year"         = "Year",
  "styear"       = "State $\\times$ Year",
  "D"            = "PESA $\\times$ Scheduled",
  "block"        = "Block",
  "blk"          = "Block"
)


# Defining vectors of estimations

ps_coef <- df_aux1$estimates[ 1:5 ]
ps_se <- df_aux1$standard.error[ 1:5 ]
cbps_coef <- df_aux1$estimates[ 6:10 ]
cbps_se <- df_aux1$standard.error[ 6:10 ]

b_ps <- c()
b_cbps <- c()
j <- 1
for ( i in 1:5){
  
  b_ps[ j ] <- ps_coef[ i ]
  b_cbps[ j ] <- cbps_coef[ i ]
  j <- j + 1
  
  b_ps[ j ] <- ps_se[ i ]
  b_cbps[ j ] <- cbps_se[ i ]
  j <- j + 1
  
}
m_table <- matrix( c( b_ps, b_cbps), 
                   ncol = 2, byrow = FALSE )
m_table <- round( m_table, 3 )

fitstat_register("n_new", function(x) summary( x )$nobs , 
                 "\\# Observations")
desc = "Dynamic Treatment Effects of PESA Adoption on 
Forest Index using Propensity Score (PS) and 
Covariate Balancing Propensity Score (CBPS)  Matching"
fn = "figure6_regs"; lab = "tab:figure6_regs";
etable( m0_ps, m1_cbps,
        style.tex = style.tex(main = "base", 
                              depvar.title = "", 
                              model.title = "", 
                              var.title = "\\midrule", 
                              slopes.title = "\\midrule \\emph{Time Trends}", 
                              yesNo = c("$\\checkmark$", "")),
        signif.code = NA,
        fixef_sizes = T, 
        fixef_sizes.simplify = F,
        headers = c( "PS Matching", "CBPS Matching" ),
        fitstat = c( "n_new" ), 
        tex = TRUE,
        dict = treatmap,
        label = lab,
        title = glue::glue("{desc}"),
        depvar = FALSE,
        extralines=list(
          "PESA X Scheduled" = m_table[ 1, ] ,
          " " = m_table[ 2, ] ,
          "Lead Year 1 X PESA X Scheduled" = m_table[ 3, ] ,
          " " = m_table[ 4, ] ,
          "Lead Year 2 X PESA X Scheduled" = m_table[ 5, ] ,
          " " = m_table[ 6, ] ,
          "Lead Year 3 X PESA X Scheduled" = m_table[ 7, ] ,
          " " = m_table[ 8, ] ,
          "Lead Year 4 X PESA X Scheduled" = m_table[ 9, ] ,
          " " = m_table[ 10, ] ,
          "\\midrule \\emph{Summary Statistics}" = c( "", "" ),
          "Mean Pre-Y (Non-Sch)"= c( rep( controls_pre1, 2) ),
          "Mean Pre-Y (Sch )"   = c( rep( treat_pre1, 2 ) ),
          "Dataset"     = c( rep( "VCF", 2 ) ),
          "Timespan"    = c( rep( "1995-2017", 2 ) )
        ),
        file = "main_figure6_tablePanelD.tex", 
        replace = TRUE
)

#------------------------------------------------------------------------------
